

def mul(a,b,c):
    print a,'=',b,'*',c,';'





def write_action():
    for i in range(1,4):
        for j in range(1,4):
            print 'Complex u%iu%i[3];' % (i,j)
            if i!=j:
                print 'Complex u%id%i[3];' % (i,j)
                print 'Complex d%iu%i[3];' % (i,j)
            print 'Complex d%id%i[3];' % (i,j)
            print 'Complex e%iu%i[3];' % (i,j)
            print 'Complex e%id%i[3];' % (i,j)
            print 'Complex u%ie%i[3];' % (i,j)
            print 'Complex d%ie%i[3];' % (i,j)
            print 'Complex b%iu%i[3];' % (i,j)
            print 'Complex b%id%i[3];' % (i,j)
            print 'Complex u%ib%i[3];' % (i,j)
            print 'Complex d%ib%i[3];' % (i,j)

    print 'Complex c_id[16];     //fermilabread_id(c_id);'
    print 'Complex c_u[16*4];    //fermilabread_u(c_u);'
    print 'Complex c_d[16*4];    //fermilabread_d(c_d);'
    print 'Complex c_e[16*4];    //fermilabread_e(c_e);'
    print 'Complex c_b[16*4];    //fermilabread_b(c_b);'
    print 'Complex c_uu[16*4*4]; //fermilabread_uu(c_uu);'
    print 'Complex c_ud[16*4*4]; //fermilabread_ud(c_ud);'
    print 'Complex c_du[16*4*4]; //fermilabread_du(c_du);'
    print 'Complex c_dd[16*4*4]; //fermilabread_dd(c_dd);'
    print 'Complex c_eu[16*4*4]; //fermilabread_eu(c_eu);'
    print 'Complex c_ed[16*4*4]; //fermilabread_ed(c_ed);'
    print 'Complex c_ue[16*4*4]; //fermilabread_ue(c_ue);'
    print 'Complex c_de[16*4*4]; //fermilabread_de(c_de);'
    print 'Complex c_bu[16*4*4]; //fermilabread_bu(c_bu);'
    print 'Complex c_bd[16*4*4]; //fermilabread_bd(c_bd);'
    print 'Complex c_ub[16*4*4]; //fermilabread_ub(c_ub);'
    print 'Complex c_db[16*4*4]; //fermilabread_db(c_db);'
    
    print 'for(a=0; a<4; a++) {\n'
    print 'forallsites(x) {\n'
    for i in range(4):
        print 'Up(x,%i)=U(x,%i)*psi_in(x+%i,a);' % (i,i,i)
    for i in range(4):
        print 'Dw(x,%i)=hermitian(U(x,%i))*psi_in(x-%i,a);' % (i,i,i)
    for i in range(1,4):
        print 'Ei(x,%i)=U.em(x,0,%i)*psi_in(x,a);' % (i-1,i)
    for i in range(1,4):
        if i==1:
            m=2
            n=3
        elif i==2:
            m=1
            n=3
        elif i==3:
            m=1
            n=2
        print 'Bi(x,%i)=U.em(x,%i,%i)*psi_in(x,a);' % (i-1,m,n)
    print '\n}\n'
    print 'Up.update();'
    print 'Dw.update();'
    print 'Ei.update();'
    print 'Bi.update();'

    print '\nforallsites(x) {\n'

    for i in range(1,4):
        for j in range(1,4):
            print 'u%iu%i=U(x,%i)*Up(x+%i,%i);' % (i,j,i,i,j)
            print 'd%id%i=hermitian(U(x,%i))*Dw(x-%i,%i);' % (i,j,i,i,j)
            if i!=j:
                print 'u%id%i=U(x,%i)*Dw(x+%i,%i);' % (i,j,i,i,j)
                print 'd%iu%i=hermitian(U(x,%i))*Up(x-%i,%i);' % (i,j,i,i,j)

            print 'u%ie%i=U(x,%i)*Ei(x+%i,%i);' % (i,j,i,i,j-1)
            print 'd%ie%i=hermitian(U(x,%i))*Ei(x-%i,%i);' % (i,j,i,i,j-1)
            print 'e%iu%i=U.em(x,0,%i)*Up(x,%i);' % (i,j,i,j)
            print 'e%id%i=U.em(x,0,%i)*Dw(x,%i);' % (i,j,i,j)
            if i==1:
                m=2
                n=3
            elif i==2:
                m=1
                n=3
            elif i==3:
                m=1
                n=2
            print 'u%ib%i=U(x,%i)* Bi(x+%i,%i);' % (i,j,i,i,j-1)
            print 'd%ib%i=hermitian(U(x,%i))*Bi(x-%i,%i);' % (i,j,i,i,j-1)
            print 'b%iu%i=U.em(x,%i,%i)*Up(x,%i);' % (i,j,m,n,j)
            print 'b%id%i=U.em(x,%i,%i)*Dw(x,%i);' % (i,i,m,n,j)


    for b in range(0,4):
        print 'k=4*a+%i;' % (b)
        print 'out=c_id[k]*psi(x,a);'

        for i in range(4):
            print 'k=16*a+4*%i+%i;' % (b,i)
            print 'if(c_u[k]!=0) out+=c_u[k]*Up(x,%i);' % (i)
            print 'if(c_d[k]!=0) out+=c_d[k]*Dw(x,%i);' % (i)
            if i>0:
                print 'if(c_e[k]!=0) out+=c_e[k]*Ei(x,%i);' % (i-1)
                print 'if(c_b[k]!=0) out+=c_b[k]*Ei(x,%i);' % (i-1)
        for i in range(1,4):
            for j in range(1,4):
                print 'k=64*a+16*%i+4*%i+%i;' % (b,i,j)
                print 'if(c_uu[k]!=0) out+=c_uu[k]*u%iu%i;' % (i,j)
                print 'if(c_dd[k]!=0) out+=c_dd[k]*d%id%i;' % (i,j)

                if i!=j:
                    print 'if(c_ud[k]!=0) out+=c_ud[k]*u%id%i;' % (i,j)
                    print 'if(c_du[k]!=0) out+=c_du[k]*d%iu%i;' % (i,j)

                print 'if(c_ue[k]!=0) out+=c_ue[k]*u%ie%i;' % (i,j)
                print 'if(c_de[k]!=0) out+=c_de[k]*d%ie%i;' % (i,j)
                print 'if(c_eu[k]!=0) out+=c_eu[k]*e%iu%i;' % (i,j)
                print 'if(c_ed[k]!=0) out+=c_ed[k]*e%id%i;' % (i,j)
                
                print 'if(c_ub[k]!=0) out+=c_ub[k]*u%ib%i;' % (i,j)
                print 'if(c_db[k]!=0) out+=c_db[k]*d%ib%i;' % (i,j)
                print 'if(c_bu[k]!=0) out+=c_bu[k]*b%iu%i;' % (i,j)
                print 'if(c_bd[k]!=0) out+=c_bd[k]*b%id%i;' % (i,j)
        print 'psi_out(x,%i)=out;' %(b)
    print '\n}'
    print '\n}'

header="""
class FermiFermilabActionNew {
 public:

  static void mul_Q(fermi_field &psi_out,
		    fermi_field &psi_in,
		    gauge_field &U,
		    coefficients &coeff, int parity=EVENODD) {
    
    mdp_lattice& lattice=psi_in.lattice()

    if(parity!=EVENODD)
      error("FermiFermilabAction::mul_Q\nparity must be EVENODD");

    if(psi_in.nspin!=4)
      error("FermiFermilabAction::mul_Q\ndoes not work for nspin!=4");
    if(psi_in.nc!=U.nc)
      error("FermiFermilabAction::mul_Q\nincompatible number of colors");

    if(coeff.has_key("sign")) sign=coeff["sign"];
    else sign=1;
    if(sign!=+1)
      error("FermiFermilabAction::mul_Q\nsign must be +1.0");

    if(coeff.has_key("kappa")) kappa_t=kappa_s=coeff["kappa"];
    if(coeff.has_key("zeta")) kappa_t=kappa_s/coeff["zeta"];

    if(coeff.has_key("kappa_s")) kappa_s=coeff["kappa_s"];
    if(coeff.has_key("kappa_t")) kappa_t=coeff["kappa_t"];


    if(kappa_t==0 || kappa_s==0) 
      error("FermiFermilabAction::mul_Q\nparameter kappa not assigned");

   
    if(coeff.has_key("r_t")) r_t=coeff["r_t"];
    else r_t=1;

    if(coeff.has_key("r_s")) r_s=coeff["r_s"];
    else r_s=1;

    if(coeff.has_key("c_{sw}")) cSW=coeff["c_{sw}"];
    else cSW=0;

    if(coeff.has_key("c_E")) c_E=coeff["c_E"];
    else c_E=1;

    if(coeff.has_key("c_B")) c_B=coeff["c_B"];
    else c_B=1;

    if(coeff.has_key("alpha_1")) alpha1=coeff["alpha_1"];
    else alpha1=0;
    if(coeff.has_key("alpha_2")) alpha2=coeff["alpha_2"];
    else alpha2=0;
    if(coeff.has_key("alpha_3")) alpha3=coeff["alpha_3"];
    else alpha3=0;
    if(coeff.has_key("alpha_4")) alpha4=coeff["alpha_4"];
    else alpha4=0;
    if(coeff.has_key("alpha_5")) alpha5=coeff["alpha_5"];
    else alpha5=0;
    if(coeff.has_key("alpha_6")) alpha6=coeff["alpha_6"];
    else alpha6=0;
    if(coeff.has_key("alpha_7")) alpha7=coeff["alpha_7"];
    else alpha7=0;
    if(coeff.has_key("alpha_8")) alpha8=coeff["alpha_8"];
    else alpha8=0;
    if(coeff.has_key("alpha_9")) alpha9=coeff["alpha_9"];
    else alpha9=0;

    site x(lattice);

    int a,k;

    NVector_field Up(lattice,4,nc);
    NVector_field Dw(lattice,4,nc);
    NVector_field Ei(lattice,3,nc);
    NVector_field Bi(lattice,3,nc);

    Matrix out(3,1);

"""

print header
write_action()
print '}\n};'
